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THE  VIRTUAL  CRACK  CLOSURE  TECHNIQUE: 

HISTORY,  APPROACH  AND  APPLICATIONS 

Ronald  Krueger* 

Abstract.  An  overview  of  the  virtual  crack  closure  technique  is  presented.  The  approach  used  is  discussed, 
the  history  summarized,  and  insight  into  its  applications  provided.  Equations  for  two-dimensional  quadrilateral 
elements  with  linear  and  quadratic  shape  functions  are  given.  Formulae  for  applying  the  technique  in  conjuction 
with  three-dimensional  solid  elements  as  well  as  plate/shell  elements  are  also  provided.  Necessary  modifications  for 
the  use  of  the  method  with  geometrically  nonlinear  finite  element  analysis  and  corrections  required  for  elements  at 
the  crack  tip  with  different  lengths  and  widths  are  discussed.  The  problems  associated  with  cracks  or  delaminations 
propagating  between  different  materials  are  mentioned  briefly,  as  well  as  a  strategy  to  minimize  these  problems. 
Due  to  an  increased  interest  in  using  a  fracture  mechanics  based  approach  to  assess  the  damage  tolerance  of 
composite  structures  in  the  design  phase  and  during  certification,  the  engineering  problems  selected  as  examples  and 
given  as  references  focus  on  the  application  of  the  technique  to  components  made  of  composite  materials. 

Key  words,  finite  element  analysis,  fracture  mechanics,  crack  closure  integral,  composite  structures, 
delamination,  interlaminar  fracture 

Subject  classincation.  Structures  and  Materials 

1.  Introduction.  One  of  the  most  common  failure  modes  for  composite  structures  is  delamination  [1-3].  The 
remote  loadings  applied  to  composite  components  are  typically  resolved  into  interlaminar  tension  and  shear  stresses 
at  discontinuities  that  create  mixed-mode  1,  II  and  III  delaminations.  To  characterize  the  onset  and  growth  of  these 
delaminations  the  use  of  fracture  mechanics  has  become  common  practice  over  the  past  two  decades  [4-6].  The  total 
strain  energy  release  rate,  Gt,  the  mode  I  component  due  to  interlaminar  tension,  G/,  the  mode  II  component  due  to 
interlaminar  sliding  shear,  G//,  and  the  mode  III  component,  G///,  due  to  interlaminar  scissoring  shear,  as  shown  in 
Figure  1,  need  to  be  calculated.  In  order  to  predict  delamination  onset  or  growth  for  two-dimensional  problems, 
these  calculated  G  components  are  compared  to  interlaminar  fracture  toughness  properties  measured  over  a  range 
from  pure  mode  I  loading  to  pure  mode  II  loading  [7-11].  A  quasi  static  mixed-mode  fracture  criterion  is  determined 
by  plotting  the  interlaminar  fracture  toughness,  Gc  ,  versus  the  mixed-mode  ratio,  Gi/Gj,  determined  from  data 
generated  using  pure  Mode  I  DCB  (Gi/Gt=0),  pure  Mode  II  4ENF  (G//G7=I),  and  mixed-mode  MMB  tests  of 
varying  ratios,  as  shown  in  Figure  2  for  IM7/8552  [12].  A  curve  fit  of  these  data  is  performed  to  determine  a 
mathematical  relationship  between  Gc  and  Gi/Gj-  [5].  Failure  is  expected  when,  for  a  given  mixed  mode  ratio  G// 
/Gr,  the  calculated  total  energy  release  rate,  Gj,  exceeds  the  interlaminar  fracture  toughness,  Gc.  Although  several 
specimens  have  also  been  suggested  for  the  measurement  of  the  mode  III  interlaminar  fracture  toughness  property 
[13-16],  an  interaction  criterion  incorporating  the  scissoring  shear,  however,  has  not  yet  been  established.  The 
virtual  crack  closure  technique  (VCCT)  [17-19]  is  widely  used  for  computing  energy  release  rates  based  on  results 


ICASE  Staff  Scientist,  Mail  Stop  132C,  NASA  Langley  Research  Center,  Hampton,  VA,  23681-2199,  (email:  rkrueger@icase.edu).  This 
research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NAS  1-97046  while  the  author  was  in 
residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 


1 


from  continuum  (2D)  and  solid  (3D)  finite  element  analyses  to  supply  the  mode  separation  required  when  using  the 
mixed-mode  fracture  criterion. 

Although  the  original  publication  on  VCCT  dates  back  a  quarter  century  [17],  the  virtual  crack  closure 
technique  has  not  yet  been  Implemented  Into  any  of  the  large  commercial  general  purpose  finite  element  codes  such 
as  MSC  NASTRAN,  ABAQUS,  ANSYS,  ASKA,  PERMAS  or  SAMCEF.  Currently  FRANC2D,  developed  by  the 
Cornell  Fracture  Group  (CFG)  at  Cornell  University,  appears  to  be  the  only  publically  available,  highly  specialized 
finite  element  code  that  uses  the  virtual  crack  closure  technique  [20,  21].  The  virtual  crack  closure  technique  has 
been  used  mainly  by  scientists  in  universities,  research  institutions  and  government  laboratories  and  is  usually 
implemented  in  their  own  specialized  codes  or  used  in  post-processing  routines  in  conjunction  with  general  purpose 
finite  element  codes.  Lately,  an  increased  interest  in  using  a  fracture  mechanics  based  approach  to  assess  the  damage 
tolerance  of  composite  structures  in  the  design  phase  and  during  certification  has  also  renewed  the  interest  in  the 
virtual  crack  closure  technique.  Efforts  are  underway  to  Incorporate  these  approaches  in  the  Composites  Material 
MIL- 17  Handbook^  . 

The  goal  of  the  current  paper  is  to  give  an  overview  of  the  virtual  crack  closure  technique,  discuss  the 
approach  used,  summarize  the  history,  and  provide  insight  into  its  application.  Equations  for  two-dimensional 
quadrilateral  elements  with  linear  and  quadratic  shape  functions  will  be  provided.  Formulae  for  applying  the 
technique  in  conjuction  with  three-dimensional  solid  elements  as  well  as  plate/shell  elements  will  also  be  given. 
Necessary  modifications  for  the  use  of  the  method  with  geometrically  nonlinear  finite  element  analysis  and 
corrections  required  for  elements  at  the  crack  tip  with  different  lengths  and  widths  will  be  discussed.  The  problems 
associated  with  cracks  or  delaminations  propagating  between  different  materials  (the  so-called  bl-material  interface) 
will  be  mentioned  briefly,  as  well  as  a  strategy  to  minimize  these  problems.  The  selected  engineering  problems 
shown  as  examples  and  given  as  references  will  focus  on  the  application  of  the  technique  related  to  composite 
materials  as  mentioned  above. 

2.  Background.  A  variety  of  methods  are  used  to  compute  the  strain  energy  release  rate  based  on  results 
obtained  from  finite  element  analysis.  Th&  finite  crack  extension  method  [22,  23]  requires  two  complete  analyses.  In 
the  model  the  crack  gets  extended  for  a  finite  length  prior  to  the  second  analysis.  The  method  provides  one  global 
total  energy  release  rate  as  global  forces  on  a  structural  level  are  multiplied  with  global  deformations  to  calculate  the 
energy  available  to  advance  the  crack.  The  virtual  crack  extension  method  [24-33]  requires  only  one  complete 
analysis  of  the  structure  to  obtain  the  deformations.  The  total  energy  release  rate  or  J-integral  is  computed  locally  at 
the  crack  front  and  the  calculation  only  involves  an  additional  computation  of  the  stiffness  matrix  of  the  elements 
affected  by  the  virtual  crack  extension.  The  method  yields  the  total  energy  release  rate  as  a  function  of  the  direction 
in  which  the  crack  was  extended  virtually,  yielding  information  on  the  most  likely  growth  direction.  Modifications 
of  the  method  have  been  suggested  in  the  literature  to  allow  the  mode  separation  for  two-dimensional  analysis  [34, 
35].  An  equivalent  domain  integral  method  which  can  be  applied  to  both  linear  and  nonlinear  problems  and 
additionally  allows  for  mode  separation  was  proposed  in  [36,  37].  The  methods  above  have  been  mentioned  here 

*  http://www.mill7.org/ 
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briefly  to  complement  the  background  information.  A  comprehensive  overview  of  different  methods  used  to 
compute  energy  release  rates  is  given  in  [38].  Alternate  approaches  to  compute  the  strain  energy  release  rate  based 
on  results  obtained  from  finite  element  analysis  have  also  been  published  recently  [39-41]. 

For  delaminations  in  laminated  composite  materials  where  the  failure  criterion  is  highly  dependend  on  the 
mixed-mode  ratio  and  propagation  occurs  in  the  laminate  plane,  the  virtual  crack  closure  technique  [17-19,  42]  has 
been  most  widely  used  for  computing  energy  release  rates  because  fracture  mode  separation  is  determined 
explicitely.  Recently  new  methods  to  compute  mixed-mode  energy  release  rates  suitable  for  the  application  with  the 
p-version  of  the  finite  element  method  have  also  been  developed  [43].  Some  modified  and  newly  developed 
formulations  allow  applications  independent  of  the  finite  element  analysis  and  suitable  for  boundary  element 
analysis  [21,  44]. 

2.1.  Crack  Closure  Method  Using  Two  Analysis  Steps.  Even  though  the  virtual  crack  closure  technique  is 
the  focus  of  this  paper  and  is  generally  mentioned  in  the  literature  it  appears  appropriate  to  include  a  related  method: 
The  crack  closure  method  or  two-step  virtual  crack  closure  technique.  The  terminology  in  the  literature  is  often 
inexact  and  this  two-step  method  is  sometimes  referred  to  as  VCCT.  It  may  be  more  appropriate  to  call  the  method 
the  crack  closure  method  because  the  crack  is  physically  extended,  or  closed,  during  two  complete  finite  element 
analyses  as  shown  in  Figure  3.  The  crack  closure  method  is  based  on  Irwin’s  crack  closure  integral  [45,  46].  The 
method  is  based  on  the  assumption  that  the  energy  Ai?  released  when  the  crack  is  extended  by  Afl  from  a  (Figure  3a) 
to  fl+Afl  (Figure  3b)  is  identical  to  the  energy  required  to  close  the  crack  between  location  I  and  i  (Figure  3a).  Index 

“1”  denotes  the  first  step  depicted  in  Figure  3a  and  index  “2”  the  second  step  as  shown  in  Figure  3b.  For  a  crack 
modeled  with  two-dimensional  four  noded  elements  as  shown  in  Figure  3  the  work  AE  required  to  close  the  crack 
along  one  element  side  can  be  calculated  as 

A£'=  -  [Xif  ■  Au2t  +  Zii  ■  Aw2i] 

where  and  are  the  shear  and  opening  forces  at  nodal  point  I  to  be  closed  (Figure  3a)  and  Au2i  and  Aw2,  are 
the  differences  in  shear  and  opening  nodal  displacements  at  node  I  as  shown  in  Figure  3b.  The  crack  closure  method 

establishes  the  original  condition  before  the  crack  was  extended.  Therefore  the  forces  required  to  close  the  crack  are 
identical  to  the  forces  acting  on  the  upper  and  lower  surfaces  of  the  closed  crack.  The  forces  X^  and  Zii  may  be 
obtained  from  a  first  finite  element  analysis  where  the  crack  is  closed  as  shown  in  (Figure  3a).  The  displacements 
Au2i  and  Aw2i  are  obtained  from  a  second  finite  element  analysis  where  the  crack  has  been  extended  to  its  full  length 

a-vAa  as  shown  in  Figure  3b. 

2.2.  The  Modified  Craek  Closure  Method.  The  modified,  or  virtual,  crack  closure  method  is  based  on  the 
same  assumptions  as  the  crack  closure  method  described  above.  Additionally,  however,  it  is  assumed  that  a  crack 
extension  of  Aa  from  a+Aa  (node  i)  to  a+2Aa  (node  k)  does  not  significantly  alter  the  state  at  the  crack  tip 
(Figure  4).  Therefore  the  displacements  behind  the  crack  tip  at  node  i  are  approximately  equal  to  the  displacements 
behind  the  original  crack  tip  at  node  1.  Further,  the  energy  AE  released  when  the  crack  is  extended  by  Aa  from  a+Aa 
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to  fl+2Aa  is  identical  to  the  energy  required  to  close  the  crack  between  location  i  and  k.  For  a  crack  modeled  with 
two-dimensional,  four-noded  elements,  as  shown  in  Figure  4,  the  work  AE  required  to  close  the  crack  along  one 
element  side  therefore  can  be  calculated  as 

A£'=^[2f,'  ■  ^{  +  Zi  ■  Awi] 

where  X,  and  Z,  are  the  shear  and  opening  forces  at  nodal  point  i  and  Aug  and  AWg  are  the  shear  and  opening 
displacements  at  node  t  as  shown  in  Figure  4.  Thus,  forces  and  displacements  required  to  calculate  the  energy  AE  to 

close  the  crack  may  be  obtained  from  one  single  finite  element  analysis.  The  details  of  calculating  the  energy  release 
rate  G=AE/AA,  where  AA  is  the  crack  surface  created,  and  the  separation  into  the  individual  mode  components  will 
be  discussed  in  the  following  section. 

3.  Equations  for  Using  the  Virtual  Crack  Closure  Technique.  In  the  following  paragraphs  equations  are 
presented  to  calculate  mixed-mode  strain  energy  release  rates  using  two-dimensional  finite  element  models  such  as 
plane  stress  or  plane  strain.  Different  approaches  are  also  discussed  for  the  cases  where  the  crack  or  delamination  is 
modeled  with  plate/shell  elements  or  with  three-dimensional  solids. 

3.1.  Formulae  for  Two-Dimensional  Analysis.  In  a  two-dimensional  finite  element  plane  stress,  or  plane 
strain  model,  the  crack  of  length  a  is  represented  as  a  one-dimensional  discontinuity  by  a  line  of  nodes  as  shown  in 
Figure  5.  Nodes  at  the  top  surface  and  the  bottom  surface  of  the  discontinuity  have  identical  coordinates,  however, 
are  not  connected  with  each  other  as  shown  in  Figure  5a.  This  lets  the  elements  connected  to  the  top  surface  of  the 
crack  deform  independently  form  those  connected  to  the  bottom  surface  and  allows  the  crack  to  open  as  shown  in 
Figure  5b.  The  crack  tip  and  the  undamaged  section,  or  the  section  where  the  crack  is  closed  and  the  structure  is  still 
intact,  is  modeled  using  single  nodes,  or  two  nodes  with  identical  coordinates  coupled  through  multi-point 
constraints  if  a  crack  propagation  analysis  is  desired.  This  is  discussed  in  detail  in  Appendix  A,  which  explains 
specific  modeling  issues. 

For  a  crack  propagation  analysis,  it  is  important  to  advance  the  crack  in  a  kinematically  compatible  way. 
Node-wise  opening/closing,  where  node  after  node  is  sequentially  released  along  the  crack,  is  possible  for  the  four- 
noded  element  as  shown  in  Figure  6a.  It  is  identical  to  element-wise  opening  in  this  case  as  the  crack  is  opened  over 
the  entire  length  of  the  element.  Node-wise  opening/closing,  however,  results  in  kinematically  incompatible 
interpenetration  for  the  eight-noded  elements  with  quadratic  shape  functions  as  shown  in  Figure  6b,  which  caused 
initial  problems  when  eight  noded  elements  were  used  in  connection  with  the  virtual  crack  closure  technique. 
Element-wise  opening  -  where  edge  and  mid-side  nodes  are  released  -  provides  a  kinematically  compatible 
condition  and  yields  reliable  results,  which  was  demonstrated  in  references  [4,  47,  48]  and  later  generalized 
expressions  to  achieve  this  were  derived  by  Raju  in  [18]. 

The  mode  I,  and  mode  II  components  of  the  strain  energy  release  rate,  G/  and  G//  are  calculated  for  four- 
noded  elements  as  shown  in  Figure  7a 
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where  Aa  is  the  length  of  the  elements  at  the  crack  front  and  X,  and  Z,  are  the  forces  at  the  crack  tip  (nodal  point  i). 
The  relative  displacements  behind  the  crack  tip  are  calculated  from  the  nodal  displacements  at  the  upper  crack  face 
Uf  and  Wf  (nodal  point  ()  and  the  nodal  displacements  and  w,*  at  the  lower  crack  face  (nodal  point  I  )  respecively. 

The  crack  surface  AA  created  is  calculated  as  AA=Aa‘\,  where  it  is  assumed  that  the  two-dimensional  model  is  of 
unit  thickness  “1”.  While  the  original  paper  by  Rybicki  and  Kanninen  is  based  on  heurisitic  arguments  [17],  Raju 
proved  the  validity  of  the  equation  [18].  He  also  showed  that  the  equations  are  applicable  if  triangular  elements, 
obtained  by  collapsing  the  rectangular  elements,  are  used  at  the  crack  tip. 

The  mode  I,  and  mode  II  components  of  the  strain  energy  release  rate,  G/,  Gu  are  calculated  for  eight-noded 
elements  as  shown  in  Figure  7b 

~  ~  2  Aa  -{ui—u^^  +  X  j  ■  {uff,  — 

where  Aa  is  the  length  of  the  elements  at  the  crack  front  as  above.  In  addition  to  the  forces  X,  and  Z,  at  the  crack  tip 
(nodal  point  i)  the  forces  X,  and  Z,  at  the  mid-side  node  in  front  of  the  crack  (nodal  point  j)  are  required.  The  relative 
sliding  and  opening  behind  the  crack  tip  are  calculated  at  nodal  points  £  and  £*  from  displacements  at  the  upper 

crack  face  Ug  and  Wg  and  the  displacements  Ug*  and  Wg*  at  the  lower  crack  face.  In  addition  to  the  relative 
displacements  at  nodal  points  I  and  £*  the  relative  displacements  at  nodal  points  m  and  m  *  are  required,  which  are 

calculated  from  displacements  at  the  upper  crack  face  and  w„g  and  the  displacements  and  w^t  at  the  lower 
crack  face  [18].  The  crack  surface  AA  created  is  calculated  as  AA=Afl‘l,  where  it  is  assumed  that  the  two- 
dimensional  model  is  of  unit  thickness  “1”.  The  equations  are  also  applicable  if  triangular  parabolic  elements, 
obtained  by  collapsing  the  parabolic  rectangular  elements,  are  used  at  the  crack  tip  [18]. 

The  total  energy  release  rate  Gj  is  calculated  from  the  individual  mode  components  as 

Gt  =  Gi+  Gu  +  G/// 

where  G///  =0for  the  two-dimensional  case  discussed. 

The  VCCT  proposed  by  Rybicki  and  Kanninen  did  not  make  any  assumptions  of  the  form  of  the  stresses  and 
displacements.  Therefore,  singularity  elements  are  not  required  at  the  crack  tip.  However,  others  have  proposed 
special  two-dimensional  crack  tip  elements  with  quarter-point  nodes  as  shown  in  Figure  8.  Based  on  the  location  of 

the  nodal  points  at  ^  =0.0,  0.25  and  1 .0,  these  quarter  point  elements  accurately  simulate  in  Ur  singularity  of  the 
stress  field  at  the  crack  tip.  Triangular  quarter-point  elements  are  obtained  by  collapsing  one  side  of  the  rectangular 
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elements,  as  shown  in  Figure  8b.  Due  to  the  fact  that  these  elements  are  not  readily  available  in  most  of  the 
commonly  used  finite  element  codes,  quarter-point  elements  are  not  discussed  further  and  the  interested  reader  is 
refered  to  the  literature  [18,  49-52] 

3.2.  Formulae  for  Three-Dimensional  Solids  and  Plate/Shell  Elements.  In  a  finite  element  model  made  of 
three-dimensional  solid  elements  (Figure  9a)  or  plate  or  shell  type  elements  (Figure  9b)  the  delamination  of  length  a 
is  represented  as  a  two-dimensional  discontinuity  by  two  surfaces.  The  additional  dimension  allows  to  calculate  the 
distribution  of  the  energy  release  rates  along  the  delamination  front  and  makes  it  possible  to  obtain  Gju,  which  is 
idential  to  zero  for  two-dimensional  models.  Nodes  at  the  top  surface  and  the  bottom  surface  have  identical 
coordinates  and  are  not  connected  with  each  other  as  explained  in  the  previous  section.  The  delamination  front  is 
represented  by  either  a  row  of  single  nodes  or  two  rows  of  nodes  with  identical  coordinates,  coupled  through  multi¬ 
point  constraints.  The  undamaged  section  where  the  delamination  is  closed  and  the  structure  is  intact  is  modeled 
using  single  nodes  or  two  nodes  with  identical  coordinates  coupled  through  multi-point  constraints  if  a  delamination 
propagation  analysis  is  desired.  This  is  discussed  in  detail  in  Appendix  A,  which  explains  specific  modeling  issues. 

3.2.1.  Formulae  for  Three-Dimensional  Solids.  For  convenience,  only  a  section  of  the  delaminated  area 
which  is  modeled  with  eight-noded  three-dimensional  solid  elements  is  illustrated  in  Figure  10.  The  mode  I, 
mode  II,  and  mode  III  components  of  the  strain  energy  release  rate,  G/,  G//,  and  Gm  are  calculated  as 

•  ^Li  ■  {wu  -w^f  ) 

Gl[=-Y^'^Li-{uu 

^  “  ''W  ) 


with  AA=Aa-b  as  shown  in  Figure  10  [53].  Here  AA  is  the  area  virtually  closed,  Aa  is  the  length  of  the  elements  at 
the  delamination  front,  and  b  is  the  width  of  the  elements.  For  better  identification  in  this  and  the  following  figures, 
columns  are  identified  by  capital  letters  and  rows  by  small  letters  as  illustrated  in  the  top  view  of  the  upper  surface 
shown  in  Figure  10b.  Hence,  Xu,  Yu  and  Zu  denote  the  forces  at  the  delamination  front  in  column  L,  row  i.  The 
corresponding  displacements  behind  the  delamination  at  the  top  face  node  row  f  are  denoted  Um,  Vki  and  w^i  and  at 

the  lower  face  node  row  I*  are  denoted  Uki*,  Vki*  and  as  shown  in  Figure  10.  All  forces  and  displacements  are 

obtained  from  the  finite  element  analysis  with  respect  to  the  global  system.  A  local  crack  tip  coordinate  system  {x 
y’,  z’),  that  defines  the  normal  and  tangential  coordinate  directions  at  the  delamination  front  in  the  deformed 
configuration,  has  been  added  to  the  illustration.  Its  use  with  respect  to  geometrically  nonlinear  analyses  will  be 
discussed  later. 

For  twenty  noded  solid  elements,  the  equations  to  calculate  the  strain  energy  release  rate  components  at  the 
element  corner  nodes  (location  Li)  as  shown  in  Figure  1 1  are 
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G/=- 


G//  = 


2M, 


2AA, 


2  ~  '^Kf  )  +  ^Li{v^U 

::  ^Ki{uKt  (“l^ 


G///=- 


2AA 


L 


M 


-ULf)  +  Xij  (  )  +  -  XMi{  -  Unjf  ) 

-  ^Li  )  +  ^Lj  [vLm  -  ^Lm  )  +  ^  ^M,(  -  W  )  ' 


where  AAL=Aa-b  as  shown  in  Figure  11  [54].  Here  X/a,  Yf^i  and  Z/f,  denote  the  forces  at  the  delamination  front  in 
column  K,  row  i.  The  relative  displacements  at  the  corresponding  column  K  are  calculated  from  the  displacements 
behind  the  delamination  at  the  lower  face  node  row  £*  as  Ugf*,  Vki*  and  Wki*  and  at  the  top  face  node  row  t,  as  Uk,  , 

Vki  and  Wk/  (Figure  11b).  Similar  definitions  are  applicable  in  column  M  for  the  forces  at  node  row  i  and 
displacements  at  node  row  (.  and  in  column  L  for  the  forces  at  node  row  i  and  j  and  displacements  at  node  row  I  and 

m  respectively.  Only  one  half  of  the  forces  at  locations  Ki  and  Mi  contribute  to  the  energy  required  to  virtually  close 
the  area  AA^.  Half  of  the  forces  at  location  Ki  contribute  to  the  closure  of  the  adjacent  area  AAj  and  half  of  the  forces 
at  location  Mi  contribute  to  the  closure  of  the  adjacent  area  AAm- 


The  equations  to  calculate  the  strain  energy  release  rate  components  at  the  mid-side  node  (location  Mi)  as 
shown  in  Figure  12  are  as  follows  [54,  55]. 


Gi 


1 


OA.,  -Wu’)+^ZLj{wijn  -^Lm-)  +  ZMi{wMf.  “  )  +  ^  -^Nf) 

Z/\A^  \_Z  Z  Z 


Gu  = 


2AAf^ 


-  Xu{uu  ~  )  +  2 


) 


^Lm  ~'^Lm 


)  +  XMi[uMl  +  Ni^ 


)  +  7;  ^Lj(vLm  ~  ^'Lm’)  + “  ^Mf  )  +  ^  ^Ni(  “ ''Nf  ) 
Z/aAj^  \_Z  Z  Z 


■  ^Nj(''Nm  ^Nm’  ) 


where  only  one  half  of  the  forces  at  locations  Li,  Lj  and  Ni,  Nj  contribute  to  the  energy  required  to  virtually  close  the 
area  AAm-  Half  of  the  forces  at  locations  Li  and  Lj  contribute  to  the  closure  of  the  adjacent  area  AA^  and  half  of  the 
forces  at  locations  Ni  and  Nj  contribute  to  the  closure  of  the  adjacent  area  AAq-  Additional  information  with  respect 
to  different  equations  for  solid  elements  are  given  in  the  literature  [19,  55-60] 


3.2.2.  Formulae  for  Plate/Shell  Elements.  The  use  of  four-noded  plate  elements  to  calculate  the  mixed¬ 
mode  strain  energy  release  was  demonstrated  in  [56].  Equations  for  the  use  of  VCCT  in  conjunction  with  four-  and 
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nine-noded  plate  elements  were  derived  by  Raju  [61].  Two  techniques  to  tie  the  nodes  at  the  delamination  front  (row 
i  in  Figues  13  and  14)  were  discussed.  For  nodes  of  the  upper  and  lower  surfaces,  which  are  joined  at  the  crack  tip, 
the  first  method  enforced  the  compatibility  of  translational  and  rotational  dregrees-of-freedom  for  nodes  with 
identical  coordinates.  In  the  second  technique,  only  the  compatibility  of  the  translational  dregrees-of-freedom  were 
enforced.  The  second  technique  allowed  the  nodes  to  have  independent  rotations.  For  the  configurations  and  loads 
considered,  the  second  technique  appeared  to  provide  preferable  constraints  yielding  accurate  strain  energy  release 
rates  [61-63].  Therefore  only  the  equations  for  this  method  are  mentioned  in  this  paper. 

For  the  four-noded  rectangular  plate  element,  as  shown  in  Figure  13,  the  equations  to  calculate  the  strain 
energy  release  rate  components  at  the  element  corner  nodes  (location  Li)  are 

■(««  ““Lf) 

where  AA=Aa(b2+b2)/2  is  the  crack  surface  closed  as  shown  in  Figure  13b  [61]. 

The  equations  to  calculate  the  strain  energy  release  rate  components  for  a  nine-noded  plate  element  as  shown 
in  Figure  14  are  for  the 


GhiIl  = 


[Yu{vLi  -VLf)  +  YLj{vLm  -^Lm')] 


GiIi\m  -  2AA^  +  “''Mm*)] 

G///|^=-2^  [^Ni{vNf.  -'’w)  +  %('’A'm  “'Wm*)] 


as  given  in  [61,  62].  Here  indices  L,  M,  N  denote  the  column  location,  as  shown  in  Figure  14b  and 


AA/_  =  —  Afl  •(^2 +^3)’  AA^  = -- Aa- (fo  1+1) 2) 

are  the  equivalent  crack  surfaces  apportioned  to  corner-  and  midside-crack  front  nodes,  respectively.  The  equivalent 
crack  surfaces  are  obtained  by  assuming  that  the  strain  energy  release  rate  components  are  constant  across  the  width 
of  an  element  [61]. 

For  eight  noded  plate  elements  the  Mm-terms  are  equal  to  zero  and  the  equations  at  column  M  are  reduced  to 

=  2AA^  [^Mi{wM( 

^Ii\m  ^  ~7AA^j  [^m(“m^  “  “Mf  )] 

=  2AA^  Me  -'’Mt)] 


with  the  equations  for  columns  L  and  N  unaltered. 

Built-up  structures  are  traditionally  modeled  and  analyzed  using  plate  or  shell  finite  elements  to  keep  the 
modeling  and  computational  effort  affordable.  Computed  mixed  mode  strain  energy  release  rate  components, 
however,  depend  on  many  variables  such  as  element  order  and  shear  deformation  assumptions,  kinematic  constraints 
in  the  neighborhood  of  the  delamination  front,  and  continuity  of  material  properties  and  section  stiffness  in  the 
vicinity  of  the  debond  when  delaminations  or  debonds  are  modeled  with  plate  or  shell  finite  elements  [61,  62].  For 
example,  in  reference  [62]  mesh  refinement  studies  showed  that  computed  G/,  G//,  and  Gui  did  not  converge  when 
the  structure  above  and  below  the  plane  of  delamination  was  modeled  with  plate  elements  with  different  section 
properties  (thickness  or  layup).  A  comparison  of  computed  mixed  mode  strain  energy  release  rates  obtained  from 
plate  models  with  values  computed  from  three-dimensional  models  showed  differences  in  results  near  the  free  edges 
of  the  structure  where  the  stress  state  is  three-dimensional  [64].  These  problems  may  be  avoided  by  using  three- 
dimensional  models.  Furthermore,  three-dimensional  analyses  are  required  when  matrix  cracks  and  multiple 
delaminations  need  to  be  modeled  at  different  ply  interfaces.  Since  many  layers  of  brick  elements  through  the 
thickness  are  often  necessary  to  model  the  individual  plies,  the  size  of  finite  element  models  required  for  accurate 
analyses  may  become  prohibitively  large.  For  future  detailed  modeling,  the  shell/3D  modeling  technique  offers  great 
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potential  for  saving  modeling  and  computational  effort  because  only  a  relatively  small  section  in  the  vicinity  of  the 
delamination  front  needs  to  modeled  with  solid  elements  [65]. 

3.3.  Formulae  for  Geometrically  Nonlinear  Analysis.  For  geometric  nonlinear  analysis  where  large 
deformations  may  occur,  both  forces  and  displacements  obtained  in  the  global  coordinate  system  need  to  be 
transformed  into  a  local  coordinate  system  (x',  z')  which  originates  at  the  crack  tip  as  shown  in  Figure  15.  The  local 
crack  tip  system  defines  the  tangential  (x',  or  mode  II)  and  normal  (z',  or  mode  I)  coordinate  directions  at  the  crack 
tip  in  the  deformed  configuration  as  shown  in  Figure  1 5b  for  the  two-dimensional  case.  The  vector  through  nodes  i 
and  k  in  the  deformed  configuration  defines  the  local  x'  direction  as  shown  in  Figure  15a.  The  local  z’  direction, 
which  defines  mode  I,  is  perpendicular  to  the  local  x’  direction  which  defines  mode  II.  Forces  at  node  row  i  and 
displacements  at  node  row  (.  need  to  be  transformed  to  the  local  x’-z’  system  at  the  tip  as  shown  in  Figure  15b.  The 

equations  to  calculate  the  mixed-mode  energy  release  rate  components  remain  the  same  as  before,  with  forces  and 
displacements  now  expressed  in  the  local  system.  For  the  two-dimensional  eight-noded  quadrilateral  element  with 
quadratic  shape  functions  this  yields 


G/  = 

^  ■\z' 

2Aa  ^  ‘ 

G//  = 

_ L  .[y 

£  —U  X  j  ‘  {Ufj2  ~  *  j  j 

2Aa  ^ 

where  X’, ,  Z’,  and  X’, ,  Z’,  are  the  forces  at  the  crack  tip  (nodal  point  i)  and  in  front  of  the  crack  (nodal  point  j)  in  the 
local  crack  tip  system.  The  relative  sliding  and  opening  behind  the  crack  tip  are  calculated  at  nodal  points  I  and 

t  from  the  transformed  displacements  at  the  upper  crack  face  u  \  and  w\  and  the  displacements  and  at  the 
lower  crack  face.  Additionally  to  the  relative  displacements  at  nodal  points  I  and  t ,  the  relative  displacements  at 

nodal  points  m  and  m  are  required,  which  are  calculated  from  displacements  at  the  upper  crack  face  u  and  w  and 
the  displacements  and  w’„t  at  the  lower  crack  face.  Three-dimensional  analysis  additionally  requires  the 
definition  of  the  tangential  (y',  or  mode  III)  coordinate  direction  which  will  be  discussed  in  section  3.5. 


3.4.  Corrections  for  Elements  with  Different  Lengths  or  Widths  at  the  Crack  Tip. 

3.4.1.  Correcting  for  Elements  with  Different  Lengths  at  the  Crack  Tip.  All  equations  in  previous 
paragraphs  have  been  derived  under  the  assumption  that  the  element  lengths  Aa  for  the  element  in  front  of  the  crack 
tip  and  behind  are  identical.  Once  automatic  mesh  generators  are  used  to  create  complex  models,  the  ideal  case  of 
identical  element  length  can  no  longer  be  assumed  and  corrections  are  required.  In  their  original  paper  Rybicki  and 

Kanninen  [17]  use  the  l/tTr  singularity  of  the  stress  field  at  the  crack  tip  to  derive  the  corrected  equations.  A  sketch 
of  a  crack  tip  modeled  with  two-dimensional  finite  elements  of  unequal  length  is  shown  in  Figure  16.  The  forces  X„ 
Z,  at  the  crack  tip  (nodal  point  i)  calculated  for  an  element  length  Aaz  are  known  from  the  finite  element  analysis. 
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Required  for  the  virtual  crack  closure  technique,  however,  are  the  forces  X,  ,  Z,  matching  the  relative  displacements 
at  node  I  behind  the  crack  tip,  which  have  been  calculated  for  an  element  length  Aa;. 

The  stress  tip  field  at  the  crack  tip  can  be  expressed  as  [46] 


1  dX  dX 

rr  dA  b- dr 


where  b  is  the  element  width  or  thickness  and  is  the  undisturbed  far  field  stress  and  o(r)  is  the  stress  in  front  of 
the  crack,  which  is  a  function  of  the  distance  r  from  the  crack  tip.  The  forces  at  the  crack  tip  for  element  length  Afl/ 
and  Aa2.  are  obtained  through  integration 

dr  1/2 

Xi=  1  ^-y-2=2 -aoo  •  Aflf' ^ 

0  r 

Xi  =  2-  b  Ooo  ■  Xa2  ^  ^ 


A  relationship  between  the  forces  can  be  derived  which  only  depends  on  the  length  of  the  elements  in  front 
and  behind  the  crack  tip 

V/2 


Aflj 

y 


The  calculation  of  the  crack  opening  force  Z,  is  done  accordingly.  With  the  relationship  of  the  forces 
established,  the  required  forces  X,,  Z,  may  be  substituted  with  the  forces  X„  Z,  obtained  from  finite  element 
analysis,  yielding  the  corrected  equations  for  the  energy  release  rate  components 

xl/2 


G/  =  - 


— —  -Z,-  -{w i  -  Wf  1 

2Aaj  ^  ’  ^Afl2  j 


xl/2 


Aa- 


A  different  approach  for  correcting  the  equations,  which  does  not  depend  on  the  assumption  of  the  l/t/r 
singularity  of  the  stress  field  at  the  crack  tip,  is  depicted  in  Figure  17  for  the  two-dimensional  case.  The  different 
element  lengths  are  accounted  for  by  correcting  the  displacements  behind  the  crack  tip  (at  node  f),  which  were 

computed  for  a  length  Afl/,  to  match  the  forces  X„  Z,  at  the  crack  tip  (nodal  point  i),  which  were  computed  for  an 
element  length  Afl/.  The  displacements  are  adjusted  by  taking  into  account  the  shape  functions  of  the  elements  or 
approximated  by  simple  linear  interpolation.  For  elements  with  linear  shape  functions,  where  the  displacement  vary 
linearly  along  the  element  edge,  both  methods  are  identical.  The  displacement  at  locations  I  and  I  are  calculated 
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using  linear  interpolation  for  Afl/  >  Aa2  as  shown  in  Figure  17a  and  linear  extrapolation  is  used  for  Afl;  <  Afl2  as 
shown  in  Figure  17b  yielding 


Afl2 

Afl, 


and 


Afl2 

Aflj 


The  calculation  of  the  crack  opening  displacement  w  is  done  accordingly.  With  the  relationship  of  the 
displacements  established,  the  required  displacements  at  locations  (.  and  (.  may  be  substituted  with  the 
displacements  f  and  t  obtained  directly  from  finite  element  analysis,  yielding  the  corrected  equations  for  the  energy 

release  rate  components 


G/  = 


1  /  \ 


2A(22  ^  ^  Aflj 


1  /  \  A^J  2 


The  method  first  described  imposes  an  analytical  relationship  based  on  the  singularity  of  the  stress 

field  at  the  crack  tip.  However,  the  second  method  is  less  restrictive  because  the  results  only  depend  on  the  finite 
element  discretisation  at  the  crack  tip. 


3.4.2.  Correcting  for  Elements  with  Different  Widths  Along  the  Delamination  Front.  For  the  derivation 
of  the  equations  in  the  previous  paragraphs  it  had  been  assumed  that  the  front  is  straight  and  that  element  widths  b 
remain  constant  along  the  front.  Meshing  of  arbitrarily  shaped  delamination  front  contours  will  however  cause 
element  length  Aa  and  width  b  to  vary  along  the  front.  Therefore,  a  two-dimensional  representation  of  a  crack  or 
delamination  in  a  plate/shell  or  a  three-dimensional  solid  finite  element  models  also  requires  a  correction  accounting 
for  differences  in  element  widths  b  along  the  front.  For  components  modeled  with  four-noded  plate/shell  type 
elements  or  eight-noded  solid  brick  elements  the  correction  is  straight  forward  as  shown  in  Figure  18.  A  variation  of 
only  the  element  length,  Aa,  yields  equations  equivalent  to  the  two-dimensional  case  discussed  earlier 


G/  = 


1  /  \  ^Cl'2 


2Aa2  •  b 

1  /  \  Ad^ 


1  /  \  Ad  2 


The  additional  variation  of  the  element  width  b  requires  the  separate  calculation  of  the  contributing  element 
surfaces  AAi  =  l/2-Afl2'^i  and  t^2=\l2- Isa^b^  yielding  the  corrected  equations  for  the  energy  release  rate  components 
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A  simple  equivalent  expression  for  twenty-noded  solid  brick  elements  based  on  the  equations  presented  in 
section  3.2.1  is  not  available.  The  set  of  equations  given  in  section  3.2.2  for  the  higher  order  plate  elements  [61]  may 
be  used  since  the  change  of  element  width  is  already  accounted  for  in  the  expressions  for  the  mixed-mode  strain 
energy  release  rates.  The  variation  In  element  length,  Afl,  may  be  compensated  as  described  in  the  section  above. 

3.5.  Procedures  for  Arbitrarily  Shaped  Delamination  Front  Contours.  The  equations  presentend  in  the 
previous  paragraphs  were  derived  with  the  assumption  that  the  delamination  front  is  straight.  For  a  straight  front  as 
shown  in  Figures  10-14  and  18  the  definition  of  the  modes  is  intuitive  and  constant  for  the  entire  front:  Mode  I  is 
caused  by  the  out  of  plane  crack  opening,  mode  II  by  the  shear  perpendicular  to  the  straight  delamination/crack  front 
and  mode  III  by  the  shear  component  tangential  to  the  front.  For  an  arbitrarily  shaped  front  the  mode  definition 
constantly  changes  along  the  contour  as  shown  in  Figure  19.  A  local  crack  tip  coordinate  system  therefore  needs  to 
be  defined  at  each  nodal  point  along  the  front  [66].  The  vector  through  nodes  i  and  k  in  the  deformed  configuration 
defines  the  local  x  direction  as  shown  for  one  sample  node  in  Figure  19.  The  secant  through  adjacent  nodes  defines 
the  y’  direction  and  mode  III.  The  local  plane  of  delamlnation  is  defined  by  and  y’ .  The  local  z  direction,  which 
defines  mode  I,  is  perpendicular  to  the  plane  of  delamination.  Finally,  the  local  x’  direction  which  defines  mode  II,  is 
perpendicular  to  the  y’  and  z  ’  plane.  Forces  at  node  row  i  and  displacements  at  node  row  i  need  to  be  transformed  to 

the  local  x’-y’-z’  system  at  the  tip.  The  transformed  forces  and  displacements  are  used  In  the  previously  derived 
equations  to  calculate  the  mixed  mode  energy  release  rates. 

For  the  procedure  described  above  it  is  assumed  that  the  mesh  remains  normal  to  the  delamination  front. 
Modern  pre-processing  software  allows  the  modeling  of  almost  any  complex  configuration,  these  programs, 
however,  were  not  developed  to  guarantee  the  normality  of  the  mesh  and  the  delamlnation  front.  The  effect  of  a  lack 
of  normality  at  the  crack  front  on  the  computed  energy  release  rates  was  studied  in  [67].  It  was  shown  that  the 
standard  formulations  for  VCCT  were  not  able  to  extract  accurate  values  from  models  that  did  not  have  normality  at 
the  crack  front  when  compared  to  reference  solutions.  It  was  also  found  that  an  Increased  variation  from  the  normal 
condition  yielded  a  greater  discrepancy.  The  formulation  of  an  extraction  method  for  VCCT  that  yields  accurate 
energy  release  rates  for  models  that  lack  normality  is  given  in  [67]. 

3.6.  Suggested  Solutions  for  Delaminations  with  Sharp  Corners.  Delaminations  with  sharp  corners  -  as 
shown  in  Figure  20a  -  pose  a  problem  when  computing  mixed  mode  energy  release  rates  as  the  separation  into  in¬ 
plane  shear  (mode  II)  and  tearing  (mode  III)  is  not  defined.  Ideally  sharp  corners  generally  do  not  exist  so  that 
modeling  a  rounded  corner,  as  shown  in  Figure  20b,  appears  to  be  an  acceptable  alternative.  The  mode  separation  is 
well  defined  at  the  rounded  corner  and  the  procedure  described  in  the  previous  paragraph  may  be  applied  to  define 
the  appropriate  local  crack  tip  coordinate  system. 
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The  method  suggested  has  been  applied  successfully  to  a  specimen  with  an  embedded,  square  delamination 
as  shown  in  Figure  21.  The  method  results  in  an  increase  of  nodes  and  elements  and  the  model  may  become  large.  A 
recently  suggested  modified  approach  uses  stair  stepped  instead  of  smoothed  fronts  and  thus  avoids  an  increase  in 
model  size  [68]. 

4.  Dealing  with  the  Problems  at  a  Bi-Material  Interfaee.  Previous  investigations  have  shown  that  care 
must  be  exercised  in  interpreting  the  values  for  G/,  Gn  and  G///  obtained  using  the  virtual  crack  closure  technique  for 
interfacial  delaminations  between  two  orthotropic  solids  [69-72].  Mathematical  solutions  of  the  near  crack  tip  field 
indicate  that  stresses  start  to  oscillate  in  the  immediate  vicinity  of  the  tip  when  crack  growth  occurs  at  interfaces 
between  materials  with  dissimilar  properties  as  shown  in  Figure  22a.  The  consequence  is  that  the  mixed  mode  ratio 
is  undefined  when  the  virtual  crack  closure  length  \a  goes  to  zero.  For  crack  growth  and  delamination  propagation 
in  composite  materials,  this  phenomenon  has  to  be  considered  as  the  delamination  is  rarely  located  at  an  interface 
between  two  plies  of  identical  orientation.  One  way  to  circumvent  this  problem  is  the  introduction  of  an  artificial 
thin  resin  rich  layer  which  is  assumed  to  exist  between  the  plies  [69,  73].  Delamination  propagation  in  this  case 
occurs  in  a  homogenious  material  and  the  above  mentioned  problem  does  not  exist.  Although  this  technique 
circumvents  the  issue,  it  requires  larger  models  with  significant  refinement  in  the  thin  resin  layer.  Several  other 
papers  have  addressed  the  issue  of  the  stress  oscillation  near  the  crack  tip  in  the  past  [74-83]. 

It  was  shown  in  the  literature  [70]  that  finite  values  for  the  virtual  crack  closure  length,  e.g.  Aa/a  >0.05  (Aa: 
virtual  crack  closure  length,  a:  crack  length),  result  in  nearly  constant  mixed  mode  ratios.  In  addition  it  was  verified 
in  [71,  72]  that  evaluating  the  energy  release  rates  via  nodal  forces  and  displacements  in  the  finite  element 
procedure,  yields  similar  results  to  the  analytical  evaluation  of  the  crack  closure  integral  based  on  the  near  tip  fields. 
For  larger  values  of  Aa  the  energy  release  rate  components  were  found  to  be  nearly  constant,  whereas  for  very  small 
Aa  they  are  functions  of  element  length  as  sketched  in  Figure  22b.  The  total  energy  release  Gt=Gj+Gi]+Gui, 
however,  converges  to  a  constant  value  as  shown  in  [62].  The  Aa/a  values  where  the  oscillation  begins  is  connected 
with  the  value  of  y  found  in  the  expression  for  the  crack  tip  singularity  which  is  given  by  It  was  also  shown 

that  evaluating  energy  release  rates  by  using  larger  values  of  Aa/a,  where  the  energy  release  rate  components 
become  stationary  yields  similar  results  to  values  obtained  by  an  analysis  performed  assuming  a  resin  rich  region 
mentioned  above  [69].  However,  the  energy  release  rates  evaluated  with  a  larger  Aa  are  insensitive  to  material 
inhomogeneities  that  exist.  This  is  a  prerequisite  for  an  energy  release  rate  that  is  to  be  used  as  fracture  criterion  in  a 
real  situation. 

For  the  virtual  crack  closure  technique,  the  energy  release  rates  are  defined  as  the  virtual  crack  closure 
integral  over  a  finite  crack  closure  length.  This  crack  closure  length  corresponds  to  the  lengths  of  the  elements 
adjacent  to  the  crack  front.  This  element  length,  Aa,  must  be  chosen  small  enough  to  assure  a  converged  FE  solution 
but  large  enough  to  avoid  oscillating  results.  The  approach  used  must  be  consistent  when  the  definition  of  the  energy 
release  rates  used  for  fracture  predictions  as  well  as  that  employed  for  material  characterisation.  This  does  not  imply 
that  the  material  tests  must  be  evaluated  by  FE-models,  but  it  should  be  established  that  the  data  reduction  scheme  is 
in  agreement  with  the  definition  of  a  finite  crack  closure  length.  Consequently,  it  had  been  suggested  to  use  element 
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lengths  at  the  crack  tip  in  such  a  manner  that  the  computed  results  are  insensitive  to  the  variation  of  the  element 
length  Aa  at  the  crack  tip  as  sketched  in  Figure  22b. 


Upper  and  lower  bounds  may  be  assumed  for  practical  applications  as  shown  in  Figure  22c.  The  element  size 
(length  and  height)  should  not  be  less  than  1/10  of  a  ply  thickness,  h,  which  corresponds  to  the  diameter  of  two 
carbon  tows  in  the  carbon/epoxy  material  modeled  as  shown  in  Figure  22c.  For  smaller  element  sizes  the  assumption 
of  modeling  each  ply  as  an  orthotropic  continuum  is  no  longer  valid.  The  ply  thickness  was  suggested  as  a  practical 
upper  limit  for  element  length  and  height  because  larger  elements  would  require  smearing  of  the  different  layer 
properties  over  one  element  as  shown  in  Figure  22c  [84].  Smeared  or  homogenized  properties  would  result  in  altered 
properties  at  the  interface  where  the  energy  release  rates  are  calculated.  Variations  in  mode  mixety  between  these 
upper  and  lower  bounds  are  typically  very  small  and  should  prove  acceptable  for  practical  applications. 

In  a  previous  investigation,  mixed  mode  energy  release  rates  were  computed  for  an  ENF  specimen  with 
multidirectional  layup,  where  the  delaminated  interface  was  located  between  a  +30°  and  -30°  ply  [85].  A  study 
indicated  that  computed  energy  release  rates,  G//  and  Gm,  do  not  exhibit  a  significant  variation  with  mesh  refinement 
as  shown  in  Table  1. 


Table  1:  Dependence  of  mixed-mode  ratio  in  bimaterial  interface  on  element  size  at  the  crack  tip 


Element  length  Aa  in 

mm 

Relative  element 

size  Aa/h 

Relative  crack 

closure  length  Aa/a 

Mode  ratio  Gu  IGt 

Mode  ratio  GuilGj 

0.0625 

0.492 

0.001969 

0.92048 

0.079509 

0.03125 

0.246 

0.000984 

0.92045 

0.079548 

0.015625 

0.123 

0.000492 

0.92020 

0.079792 

0.078125 

0.0615 

0.000246 

0.91574 

0.084262 

0.00390625 

0.0307 

0.000123 

0.91562 

0.084377 

Layup  [+30/-30/30/-30/30/30/-30/30/-30/-30/30/t-30/30/30/-30/30/-30/-30/30/-30/30/+30],C12K/R6376tape, 


delamination  length  a=31.75  mm,  ply  thickness  /t=0.127  mm 


In  another  study  the  influence  of  the  mesh  size  at  the  delamination  tip  was  investigated  for  a  Single  Leg 
Bending  (SLB)  specimen  with  multidirectional  layup,  where  the  delaminated  interface  was  located  between  a  +30° 
and  -30°  ply  [86].  The  three-dimensional  model  of  the  SLB  specimen  is  shown  in  Figure  23.  Along  the  length  of  the 
model,  a  refined  mesh  of  length,  c,  was  used  in  the  vicinity  of  the  delamination  front.  The  influence  of  mesh  size  on 
computed  mixed  mode  strain  energy  release  rates  was  studied  by  keeping  the  length  of  the  refined  zone,  c,  constant, 
and  increasing  the  number  of  elements,  n,  in  this  zone  as  shown  in  the  detail  of  Figure  23.  The  corresponding  values 
of  relative  element  size,  Aa/h,  and  relative  crack  closure  length,  Aa/a,  are  given  in  Table  2.  The  influence  of  mesh 
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refinement  on  the  mode  I  strain  energy  release  rate  distribution  across  the  width  is  moderate  and  only  very  long 
elements  (n=3,  Aa=c/n=l  mm)  need  to  be  avoided  as  shown  in  Figure  24.  This  is  confirmed  by  the  mode  II  and 
mode  III  distributions  as  shown  in  Figures  25  and  26  where  the  mode  II  strain  energy  release  rate  is  fairly  constant 
across  almost  the  entire  width  of  the  specimen  and  peaks  near  the  edges  accompanied  by  local  mode  III  contribution. 
The  distribution  of  the  mixed  mode  ratio  Gi/Gn  is  shown  in  Figure  27.  For  the  range  studied  (n=3  up  to  48),  there  is 
only  a  small  dependence  of  computed  mixed  mode  ratio  on  element  size  \a. 


Table  2:  Mesh  size  at  the  crack  tip  for  SLB  specimen  with  multidirectional  layup. 


length  c  of  refined 

section  in  mm 

Number  of  elements 

n 

Element  length  i^a  in 

mm 

Relative  element 

size  Aa/h 

Relative  crack 

closure  length  Aa/a 

3.0 

3 

1.0 

7.874 

0.029 

3.0 

6 

0.5 

3.937 

0.0145 

3.0 

12 

0.25 

1.9685 

0.0073 

3.0 

24 

0.125 

0.984 

0.00365 

3.0 

48 

0.0625 

0.492 

0.00183 

Layup  [+30/0/-30/0/30/04/30/0/-30/0/-30/30/T-30/30/30/0/30/0/-30/04/-30/0/30/0/+30],C12K/R6376tape 
delamlnatlon  length  <2=34.3  mm,  ply  thickness  /j=0.127  mm 


The  results  from  the  above  studies  of  the  multidirectional  ENF  and  SLB  specimens  confirm  the  suggestion 
made  earlier  for  lower  bounds  {Aa/h  >0.1)  and  upper  bounds  (Aa/h  <1.0)  of  element  lengths  to  be  used. 


5.  Application  of  the  Virtual  Crack  Closure  Technique  to  Engineering  Problems.  Selected  engineering 
problems  where  the  virtual  crack  closure  technique  is  applied  to  fracture  of  composite  structures  are  referenced  in 
this  section.  Two  examples  from  recent  studies  showing  the  application  of  VCCT  to  two-dimensional  and  three- 
dimensional  finite  element  analyses  are  discussed  in  more  detail. 

5.1.  Applying  the  Virtual  Crack  Closure  Technique  with  Two-Dimensional  Finite-Element  Analysis. 

Two-dimensional  finite  element  models  where  the  crack  or  delamination  is  simulated  as  a  line  have  been  used 
extensively  in  the  past.  In  general  two-dimensional  models  are  preferred  by  industry  due  to  the  fact  that  modeling 
time,  as  well  as  computational  time,  remains  affordable,  especially  if  many  different  configurations  have  to  be 
analyzed  during  the  initial  design  phase.  However,  the  geometry,  boundary  conditions,  and  other  properties  across 
the  entire  width  are  are  inherently  constant  in  a  two-dimensional  finite  element  model.  Additionally,  a  two- 
dimensional  plane-stress  model  in  the  x-y  plane  imposes  the  out  of  plane  stresses  to  be  zero  (  =t:xz  yz  —  0) 
and  allows  the  displacement  to  be  the  free  parameter.  A  plane-strain  model  in  the  x-y  plane,  on  the  other  hand. 
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imposes  the  out  of  plane  strains  to  be  zero  =Yjcz  yz  =0)’  which  excessively  constrains  the  plies.  The  effect 
of  two-dimensional  modeling  assumptions  is  most  marked  for  45°  plies  because  of  their  high  in-plane  Poisson’s 
ratio,  while  it  is  small  for  0°  and  90°.  The  influence  of  two-dimensional  finite  element  modeling  assumptions  on 
computed  energy  release  rates  was  studied  in  detail  in  [87].  Based  on  the  results  of  this  investigation  it  is 
recommended  to  use  results  from  plane  stress  and  plane  strain  models  as  upper  and  lower  bounds.  Two-dimensional 
models  may  also  be  used  to  qualitatively  evaluate  the  variation  of  energy  release  rates  and  mixed  mode  ratios  with 
delamination  length.  For  more  accurate  predictions,  however,  a  three-dimensional  analysis  is  required. 

Two-dimensional  finite  element  models  have  been  used  to  study  the  behavior  of  specimens  used  and 
suggested  for  fracture  toughness  testing  in  [88-95]  and  three-point  bending  specimen  [96].  losipescu  specimens  were 
studied  in  [97,  98].  The  behavior  of  edge  delaminations  was  investigated  in  [4,  99]  and  the  failure  of  composite  hat 
stringer  pull-of  specimens  was  examined  in  [100,  101].  The  failure  of  a  lap  joints  was  investigated  in  [102,  103]  and 
the  durability  of  bonded  joints  in  [104].  The  virtual  crack  closure  technique  was  also  used  to  investigate  the 
facesheet  delaminating  from  the  sandwich  core  [105-108]  as  well  as  delamination  buckling  [109,  110]. 
Delamination  initiation  from  ply  drops  in  general  was  studied  in  [111]  while  the  initiation  specific  to  rotorcraft 
flexbeams  was  investigated  in  [1 12,  113]. 

To  illustrate  the  application  of  VCCT  to  structural  delamination,  an  example  is  taken  from  the  study  of 
skin/stringer  debond  failure  as  shown  in  Figure  28  [114].  The  specimens  consisted  of  a  tapered  laminate, 
representing  the  stringer,  bonded  onto  a  skin  as  shown  in  Figure  28a.  An  IM7/8552  graphite/epoxy  system  was  used 
for  both  the  skin  and  flange.  The  skin  was  made  of  prepreg  tape  and  had  a  nominal  ply  thickness  of  0. 142  mm  and  a 
[45/-45/0/-45/45/90/90/-45/45/0/45/-45]  lay-up.  The  flange  was  made  of  a  plain-weave  fabric  of  0.208  mm  nominal 
ply  thickness.  The  flange  lay-up  was  [45/0/45/0/45/0/45/0/45]f,  where  the  subscript  “f’  denotes  fabric,  “0” 
represents  a  0°-90°  fabric  ply  and  “45”  represents  a  0°-90°  fabric  ply  rotated  by  45°.  Quasi-static  tension  tests  were 
performed  in  a  servohydraulic  load  frame  where  the  specimens  were  mounted  in  hydraulic  grips  with  a  gage  length 
of  101.6  mm  as  shown  in  Figure  28b.  The  value  of  the  damage  onset  load,  P  was  averaged  from  four  tests  and 
determined  to  be  17.8  kN  with  a  coefficient  of  variation  of  8.9%.  The  tests  were  terminated  when  the  flange 
debonded  from  the  skin.  Micrographic  investigations  showed  that  at  corners  2  and  3,  a  delamination  formed  at  the 
top  45°/-45°  skin  ply  interface.  The  initial  crack  was  modeled  at  the  location  suggested  by  the  microscopic 
investigation  as  shown  in  the  detail  in  Figure  28b.  Finite  element  solutions  were  obtained  using  the  commercial 
ABAQUS®/Standard  finite  element  software.  Eight-noded  quadrilateral  plane-stress  (CPS8R)  elements  with 
quadratic  shape  functions  and  a  reduced  (2x2)  integration  scheme  were  utilized  for  the  geometric  nonlinear  analyses. 
For  this  investigation,  the  delamination  growing  between  the  skin  top  45°  and  -45°  plies  was  extended  by  releasing 
multi-point  constraints  at  the  crack  tip  and  in  front  of  the  crack  tip.  During  a  series  of  nonlinear  finite  element 
analyses,  strain  energy  release  rates  were  computed  at  each  tip  location  for  the  loads  applied  in  the  experiments.  A 
critical  energy  release  rate,  Gc,  needs  to  be  determined  to  predict  delamination  onset.  This  critical  G  is  generally 
identified  based  on  the  shape  of  the  total  energy  release  rate  versus  delamination  length  curve,  which  is  determined 
through  analysis  as  shown  in  Figure  29.  The  Gj  versus  .r  curve  reached  a  maximum  (as  marked  in  Figure  29)  at 
some  virtual  delamination  length  and  then  decreased.  The  delamination  was  extended  to  a  total  simulated  length  of 
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2.2  mm  to  ascertain  that  the  peak  value  had  been  captured.  The  total  energy  release  rates  and  the  corresponding 
mixed  mode  ratios  Gj/Gj  are  plotted  in  Figures  29  and  30.  The  maximum  total  energy  release  rate,  Gj,  as  marked  in 
the  graph,  was  chosen  as  the  critical  value  to  be  used  for  the  fatigue  life  prediction  [114].  This  peak  in  Gj  also 
corresponded  to  the  maximum  mode  I  percentage,  i.e.,  minimum  value  of  Gi/Gj  in  Figure  30.  Similar  investigations 
of  skin/flange  debonding  are  reported  in  [115-118]. 

5.2.  Applying  the  Virtual  Crack  Closure  Technique  to  Analysis  with  Plate  and  Shell  Elements.  In  a 

finite  element  model  made  of  plate  or  shell  type  elements  the  delamination  is  represented  as  a  two-dimensional 
discontinuity  by  two  surfaces  with  identical  nodal  point  coordinates.  Flanagan  developed  a  code  based  on 
sublaminate  analysis  and  demonstrated  its  application  to  simple  fracture  toughness  specimens,  the  curved-beam  test 
and  edge  delamination  [119].  The  problem  of  a  composite  stringer  separating  from  the  skin  has  been  investigated 
using  plate  and  shell  models  in  [61,  120,  121].  Fractures  mechnics  analyses  performed  by  Glaessgen  et.al.  focussed 
on  the  debonding  of  stitched  composite  structures  [122-124].  Delamination  buckling  has  been  investigated  as  well 
using  plate/ shell  models  [56,  125-129]. 

5.3.  Applying  the  Virtual  Crack  Closure  Technique  to  Analysis  with  Solid  Finite  Elements.  Three- 
dimensional  finite  element  models  have  been  used  to  study  the  behavior  of  specimens  traditionally  used  in  fracture 
toughness  testing  [16,  23,  55,  64,  65,  85,  130-136],  three-point  bending  specimen  [137]  as  well  as  the  behavior  of 
edge  delaminations  [138-140].  Skin-stiffener  debonding  was  analysed  in  [54,  87].  Damage  in  titanium-grahite 
hybrid  lamiantes  was  investigated  in  [141]  and  Ireman  et.al.  studied  the  damage  propagation  in  composite  structural 
elements  [142].  Three-dimensional  models  were  also  used  to  study  the  growth  of  near-surface  delaminations  in 
composite  laminates  under  fatigue  loading  [143-145].  Delamination  buckling  has  been  investigated  extensively 
using  three-dimensional  models  [56,  66,  111,  146-156].  A  comparison  of  two-dimensional  and  three-dimensional 
analysis  for  skin/stringer  debond  failure  was  performed  in  [87,  157,  158]. 

The  example  chosen  is  taken  from  an  investigation  of  delamination  buckling  and  growth  as  shown  in 
Figure  31a  [143,  144,  159].  The  specimens  were  made  of  T300/914C  tape  material  with  a  [+5/+45/+5/-45/0/+85/0/- 
45/-+5/+45/-+5]  layup.  An  embedded  circular  delamination  of  10  mm  diameter  at  interface  2/3  was  assumed  to 
simulate  an  impact  damage  near  the  surface  as  shown  in  Figure  31b.  The  specimens  were  subjected  to  tension- 
compression  (R=-l)  fatigue  loading.  Stress  maxima  in  the  range  of  220-240  N/mm^  had  shown  to  yield  stable 
delamination  growth.  During  the  experiment,  the  out-of-plane  (i.e.  buckling)  deformation  was  monitored  by  Moire 
technique.  Using  numerical  post-processing  procedures,  the  size  and  shape  of  the  delaminated  sublaminate  was 
determined  from  this  information,  yielding  the  delamination  contours,  which  were  used  as  input  to  the  numerical 
model.  Additional  X-ray  photographies  (as  shown  in  Figure  31b)  were  used  to  verify  this  method.  Finite  element 
solutions  were  obtained  using  the  NOVA  finite  element  software.  Due  to  the  extensive  computation  times  seen  for 
three-dimensional  models,  a  special  layered  element  with  eight  nodes,  formulated  according  to  a  continuum  based 
three  dimensional  shell  theory,  was  used  for  the  nonlinear  simulations  [160].  Interpenetration  of  the  layers  in  the 
delaminated  region  was  prevented  by  using  a  contact  processor  that  utilizes  a  contactor  target  concept  applying  the 
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penalty  method  [161].  Figure  31a  shows  the  specimen  outline,  details  of  the  section  modeled  and  the  deformed  FE- 
model.  The  distributions  of  the  mixed  mode  energy  release  rates  along  the  front  S2  (i.e.  after  200,000  load  cycles)  as 
shown  in  Figure  31b,  are  plotted  in  Figure  32.  The  global  maximum  of  the  total  energy  release  rate  G^is  now  found 
perpendicular  to  the  loading  direction  (i=0.25  and  i=0.75)  with  a  second  maximum  occuring  in  the  loading 
direction.  The  values  for  interlaminar  shear  failure  also  reach  their  maximum  perpendicular  to  the  load  direction 
(i=0.22  and  i=0.72)  and  a  second  maximum  for  i=0  and  i=0.5.  This  also  holds  for  the  G/  distribution  where  again 
the  maximum  is  reached  perpendicular  to  the  loading  direction,  in  loading  direction  however  it  ceases  to  zero. 
Details  of  the  entire  investigation  are  discussed  in  [162]. 

6.  Concluding  Remarks.  The  increased  interest  in  using  a  fracture  mechanics  based  approach  to  assess  the 
damage  tolerance  of  composite  structures  in  the  design  phase  and  during  certification  has  renewed  the  interest  in  the 
virtual  crack  closure  technique.  Efforts  are  underway  to  incorporate  these  approaches  in  the  Composites  Material 
MIL-17  Handbook,  which  has  been  the  motivation  for  the  overview  presented. 

The  approach  used  in  the  virtual  crack  closure  technique  was  discussed.  Equations  for  the  use  of  the 
technique  in  conjunction  with  two-dimensional  elements,  three-elements  solids  as  well  as  plate/shell  elements  were 
given  and  insight  into  the  application  to  engineering  problems  was  provided.  The  paper,  however,  is  not  a 
comprehensive  literature  survey,  but  more  a  summary  and  review  of  issues  relevant  to  the  successful  application  of 
the  virtual  crack  closure  technique.  Nevertheless,  an  effort  was  made  to  pay  proper  tribute  to  the  relevant 
contributions  made  during  a  quarter  century  since  the  original  publication. 
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Appendix  A. 

Implementation  of  the  Virtual  Crack  Closure  Technique.  Currently,  the  large  commercial  finite  element 
codes  do  not  offer  the  choice  for  calculating  the  mixed  mode  energy  release  rates  using  the  virtual  crack  closure 
technique  (VCCT)  as  described  in  the  main  text.  Therefore,  the  energy  release  rate  components  G/ ,  G//  and  Gui  need 
to  be  computed  by  user  written  subroutines.  These  subroutines  may  either  interface  directly  with  the  finite  element 
code  during  its  execution,  provided  this  option  has  been  made  available  for  this  particular  software,  or  operate  as  an 
entirely  separate  post-processing  step.  In  both  cases  the  strategy  is  similar  and  attention  has  to  be  given  to  the 
modeling  of  the  discontinuity  as  shown  in  Figure  A1  and  the  calculation  of  the  mixed  mode  energy  release  rate  as 
outlined  in  the  flowchart  of  Figure  A2. 

As  mentioned  in  section  3.1,  the  crack  of  length,  a,  in  a  two-dimensional  finite  element  model  is  represented 
as  a  one-dimensional  discontinuity  by  a  line  of  nodes  as  shown  in  Figures  5  and  33.  Nodes  at  the  top  surface  and  the 
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bottom  surface  have  identical  coordinates,  however,  are  not  connected  with  each  other  as  shown  in  Figure  33a.  The 
crack  tip  is  represented  by  either  a  single  node  (Figure  33a)  or  two  nodes  with  identical  coordinates  coupled  through 
multi-point  constraints  (Figure  33b).  The  undamaged  section  or  the  section  where  the  crack  is  closed  and  the 
structure  is  still  Intact  is  modeled  using  single  nodes  (Figure  33a)  or  two  nodes  with  identical  coordinates  coupled 
through  multi-point  constraints  (Figure  33b).  It  is  generally  up  to  the  user  to  decide  how  to  model  the  crack  tip  and 
the  intact  region.  However,  the  use  of  multi-point  constraints  may  be  preferred  if  a  crack  propagation  analysis  is  to 
be  performed.  The  use  of  multi-point  constraints  may  also  be  preferable  if  the  finite  element  code  used  provides 
“forces  at  constraints”  as  a  standard  output.  This  feature  would  render  obsolete  the  step  where  element  forces  at 
nodes  have  to  be  retrieved  and  summed  to  obtain  forces  at  the  crack  tip.  The  same  logic  applies  to  a  finite  element 
model  made  of  plate/shell  type  elements  or  three-dimensional  solid  elements,  where  the  delamlnation  of  length,  a,  is 
represented  as  a  two-dimensional  discontinuity  by  two  surfaces. 

The  application  of  the  virtual  crack  closure  technique  based  on  results  from  finite  element  analysis  requires 
access  to  the  element  forces  at  nodes,  the  nodal  point  displacements  and  the  nodal  point  coordinates.  The  flow  chart, 
depicted  in  Figure  34  as  an  example,  shows  an  Independent  post-processing  procedure  where  input  data  for  the 
VCCT  was  extracted  directly  from  an  ABAQUS®  binary  result  file.  The  strategy  to  use  VCCT,  however,  is  the  same 
if  data  is  retrived  from  a  data  file  in  ASCII  format  or  if  a  user  written  subroutine  interfaces  directly  with  the  finite 
element  software  during  execution: 

1.  Establish  interface  with  finite  element  sortware  if  required. 

2.  Provide  input  interface  to  read  variable  problem  specific  external  data  such  as  crack  length.  Young’s 
modulus,  etc.,  and  control  parameters  such  as  element  type,  etc.,  which  is  generally  not  hard  coded 
into  the  subroutine. 

3.  Create  output  file  to  store  echo  print  of  control  parameters.  Input  data,  and  retrieved  data  as  well  as 
intermediate  and  final  results. 

4.  Read  element  forces,  or  forces  at  constraints,  nodal  point  coordinates,  and  nodal  displacements  from 
binary  result  file,  or  output  file  in  ASCII  format  or  data  bank  if  data  is  stored  externally. 
Alternatively  retrieve  required  data  from  specific  common  blocks  or  arrays  of  the  finite  element 
software  If  a  user  Interface  is  provided. 

5.  Store  the  data  retrieved  in  step  4  in  external  files,  data  base,  or  internal  common  blocks  and  arrays 
for  further  access  during  the  calculation. 

6.  Calculate  area  virtually  closed,  AA,  using  nodal  point  coordinates. 

7.  Calculate  local  coordinate  system  x\y’  and  z’  at  crack  tip  for  geometrically  nonlinear  analysis  and 
arbitrarily  shaped  delamlnation  front  as  discussed  in  sections  3.3,  3.5,  and  3.6. 

8. 

a.  Obtain  the  forces  at  the  crack  tip  and  in  front  of  the  crack  tip  from  the  forces  at  element 
nodes  by  summing  the  forces  at  common  nodes  from  elements  belonging  either  to  the 
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upper  or  the  lower  surface.  This  step  may  be  skipped  if  forces  at  contraints  are  available 
directly  from  the  finite  element  software. 


b.  Obtain  relative  displacements  between  the  surfaces  from  the  nodal  point  displacements. 

c.  Transform  forces  and  displacements  to  local  crack  tip  coordinate  system. 

9.  Calculate  mixed  mode  strain  energy  release  rates  using  the  equations  given  for  the  individual 
element  type  and  correct  for  different  element  lengths  and  widths  if  necessary. 

10.  Repeat  steps  6  through  9  for  all  nodes  along  the  delamination  front  and  write  results  to  output  file. 
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FIGURE  1 .  Fracture  Modes. 
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a  Aa  Aa 


(a).  First  Step  -  Crack  closed 


a  Aa  Aa 
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a  Aa  Aa 


FIGURE  4.  Modified  Crack  Closure  Method  (One  step  VCCT). 
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crack  closed 
undamaged  structure 


FIGURE  5.  Crack  modeled  as  one-dimensional  discontinuity. 
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(a).  Node  wise  crack  opening  for  four-noded  element 


a 


2nd  configuration:  incompatible  node-wise  2nd  configuration:  compatible  element-wise 

crack  opening/closure  crack  opening/closure 


(b).  Crack  opening  for  eight-noded  element 
FIGURE  6.  Kinematic  compatible  crack  opening/closure. 
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a  Aa  Aa 


Gii  =  -  Xi  (  u^-  U/  )  /  (  2Aa  ) 

(a).  Virtual  Crack  Closure  Technique  for  four-noded  element. 


a  Aa  Aa 


G|  =  -[  Zi  ( vVf-  w^* )  +  Z'j  (  Wm-  Wm* )  ]  /  (  2Aa  ) 

Gii  =  -[  Xi  (  U;  -  Ur  )  +  X'j  (  Um  -  Um* )  ]  /  (  2Aa  ) 

(b).  Virtual  Crack  Closure  Technique  for  eight-noded  element. 

FIGURE  7.  Virtual  Crack  Closure  Technique  for  two-dimensional  solid  elements. 
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(a).  Quadtrilateral  elements  with  quarter  point  nodes. 


a  Aa  Aa 


(b).  Collapsed  quarter  point  element. 

FIGURE  8.  Singularity  elements  with  quarter  point  nodes  at  crack  tip. 
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nodal  points  at 


(a).  Delamination  modeled  with  bilinear  three-dimensional  solid  elements 


delamination  front 
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(b).  Top  view  of  upper  surface  (lower  surface  terms  are  omitted  for  clarity) 

FIGURE  10.  Virtual  Crack  Closure  Technique  for  four  noded  plate/shell  and  eight  noded  solid  elements. 
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(b).  Top  view  of  upper  surface  (lower  surface  terms  are  omitted  for  clarity) 

FIGURE  1 1 .  Virtual  Crack  Closure  Technique  for  corner  nodes  in  eight  noded  plate/shell  and  twenty  noded  solid  elements. 
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local  system 


z',w',Z' 


FIGURE  12.  Virtual  Crack  Closure  Technique  for  midside  nodes  in  eight  noded  plate/shell  and  twenty  noded  solid  elements. 
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local  crack  tip 
system 


z',w',Z' 


(b).  Top  view  of  upper  surface  (lower  surface  terms  are  omitted  for  clarity) 


FIGURE  13.  Virtual  Crack  Closure  Technique  for  four  noded  plate/shell  elements. 
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(b).  Top  view  of  upper  surface  (lower  surface  terms  are  omitted  for  clarity) 


FIGURE  14.  Virtual  Crack  Closure  Technique  for  nine  noded  plate/shell  elements. 
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(b).  Forces  and  displacements  in  local  coordinate  system 
FIGURE  15:  Virtual  Crack  Closure  Technique  for  geometrically  nonlinear  analysis. 


43 


Aai  Aa2 


FIGURE  16.  Correction  for  elements  with  different  lengths  in  front  and  behind  the  crack  tip  . 
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a  Aa-i  Aa2 


FIGURE  17.  Correction  for  elements  with  different  lengths  in  front  and  behind  the  crack  tip  . 


45 


delamination  front 


delaminated 

area 


bi 


b2 


bi/2  j 

b2/2  i 


i  y'.v'.Y' 


VL^ 


ULi 


AAi 


^aa2; 


Yu 


Xu 


il - 


J » - ( 


Aai 


i  k 


Aa2 


intact  area 


x',u',X' 
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Top  view  of  upper  surface 
(lower  surface  terms  are  omitted  for  clarity) 

FIGURE  18.  Virtual  Crack  Closure  Technique  for  eight  noded  solid  elements  with  different  widths. 
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In-plane  shear 
Mode  II 


In-plane  shear 
Mode  II 


Tearing 
Mode  III 


(b).  Well  defined  mode  separation  at  rounded  corner. 
FIGURE  20.  Definition  of  mode  separation  at  sharp  corners. 
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FIGURE  21.  Detail  of  a  finite  element  mesh  with  modeled  square  insert. 
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Aa:  element  length 
h:  ply  thickness 


(b).  Dependence  of  computed  energy  release  rate  on  element  size  at  crack  tip 


1.0  h 


individual  plies 


(c).  Upper  and  lower  bounds  for  element  size  at  crack  tip 


FIGURE  22.  Bi-Material  Interface. 
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- outline  of  undeformed  configuration 

Wc:  deflection  in  the  center  of  the  specimen 
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FIGURE  23.  Three-dimensional  finite  element  model  ofSLB  specimen. 


51 


0.030 


0.025 


0.020  - 


Energy 
Release 
RateG,,  0.015 


kJ/m" 


0.010  - 


0.005  - 


0.000 


-0.4 


o 

n=3 

□ 

n=6 

o 

n=12 

A 

n=24 

V 

n=48 

o 


o 


o 


-0.2 


0 

y/B 


0.2 


o 


0.4 


FIGURE  24.  Influence  of  number  of  elements  in  refined  section  on  computed  mode  I  strain  energy 
release  rate  distribution  across  the  width  of  a  SLB  specimen. 
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FIGURE  25.  Influence  of  number  of  elements  in  refined  section  on  computed  mode  II  strain  energy 
release  rate  distribution  across  the  width  of  a  SLB  specimen. 
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FIGURE  26.  Influence  of  number  of  elements  in  refined  section  on  computed  mode  III  strain  energy 
release  rate  distribution  across  the  width  of  a  SLB  specimen. 
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FIGURE  27.  Influence  of  number  of  elements  in  refined  section  on  mode  ratio  distribution  across  the  width  of  a  SLB  specimen. 
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(f)).  Finite  element  model  of  skin/flange  specimen 
FIGURE  2% .  Skin/ flange  debonding  specimen. 
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X  coordinate,  mm 


FIGURE  29.  Computed  total  strain  energy  release  rate  for  delamination  between  top  45°/-45°  skin  plies. 


X  coordinate,  mm 

FIGURE  30.  Computed  mixed  mode  ratio  for  delamination  between  top  45°/-45°  skin  plies. 
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87.0 


(a).  Delamination  buckling  specimen  and  section  modeled  with  three-dimensional  solid  elements 


(b).  X-Ray  showing  the  initial  circular  delamination  and  detected  growth  after  200,000  load  cycles 


FIGURE  31.  Delamination  buckling  specimen. 
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FIGURE  32.  Computed  mixed  mode  strain  energy  release  rate  along  detected  delamination  front  after  200,000  load  cycles. 
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crack  closed 
undamaged  structure 


(b).  Intact  region  in  front  of  the  crack  modeled  with  double  nodes  tied  by  multi-point  contsraints 
FIGURE  Al.  Crack  modeled  as  one-dimensional  discontinuity 
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FIGURE  A2.  Flow  chart  of  routine  extractfto  calculate  strain  energy  release  rates. 
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